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Abstract 



We use molecular dynamics simulations in two dimensions to investigate 
the possibility that a core-softened potential can reproduce static and dynamic 
anomalies found experimentally in liquid water: (i) the increase in specific vol- 
ume upon cooling, (ii) the increase in isothermal compressibility upon cooling, 
and (iii) the increase in the diffusion coefficient with pressure. We relate these 
anomalies to the shape of the potential. We obtain the phase diagram of the 
system and identify two solid phases: a square crystal (high density phase), 
and a triangular crystal (low density phase). We also discuss the relation 
between the anomalies observed and the polymorphism of the solid. Finally, 
we compare the phase diagram of our model system with experimental data, 
noting especially the line of temperatures of maximum density, line of pres- 
sures of maximum diffusion constant, and line of temperatures of minimum 
isothermal compressibility. 

PACS numbers: 61.20.Gy, 61.25.Em, 65.70.+y, 64.70.Ja 
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I. INTRODUCTION 



Most liquids contract upon cooling and become more viscous with pressure. This is not 
the case for the most important liquid on earth, water. For at least 300 years it has been 
known that the specific volume of water at ambient pressure starts to increase when cooled 
below T = 4°C ||T|. It is perhaps less known that the viscosity of water decreases upon 
increasing pressure in a certain range of temperatures 0. Moreover, in a certain range 
of pressures water exhibits an anomalous increase of compressibility, and hence of density 
fluctuations, upon cooling. These anomalies are not restricted to water but are also present 
in other liquids ||-[5| . 

In order to investigate these anomalies, we utilize computer simulation of a class of 
potentials called "core softened" potentials, first introduced by Stell, Hemmer, and their 
coworkers ||. We define a core-softened potential as a spherically symmetric potential that 
has a region of negative curvature in its repulsive core 0. An example of a discrete and 
of a smooth core-softened potential is shown in Fig. |l[ Debenedetti et al. noted that a 
"softened core" can lead to a density anomaly ||, i.e., one of the anomalies found in water. 
Furthermore, ab initio calculation || and inversion of the experimental oxygen-oxygen radial 
distribution function reveals that a "core-softened" potential can be considered a realistic 
first-order approximation for the interaction between water molecules 1TD . 



Although directional bonding is certainly a fundamental issue in obtaining quantitative 
predictions for network-forming liquids like water, it could be the case that core-softened 
potentials can be the simplest framework to understand the physics of those anomalies. 
Here we demonstrate, by means of numerical simulations for d = 2, that the core-softened 
potential can lead to anomalies in the density, in the compressibility and in the viscosity. 
We also offer an explanation for the occurrence of these three anomalies in terms of the 
shape of the potential. 

The paper is organized as follows: in Sec. we define the potentials studied and in 
Sec. [Til] we describe the methods of simulations employed. In Sec. [TV] we present the results 
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about the density anomaly. In Sec. |Vj we discuss the relation of the structures in the solid 
phase and the density anomalies. In Sec. [VI] we present the results for the diffusion anomaly 
and give an explanation of such anomaly in terms of free volume. In Sec. |V11| we present 
the results on the compressibility anomaly. Finally, we present the overall phase diagram in 
Sec. IV 1111 and our conclusions and comments in Sec. ITX. 



II. DISCRETE AND SMOOTH MODELS 

A. Discrete Potential 

The core-softened potentials that we study are shown in Fig. [I] as a function of particle 
pair distance r. The discrete potential is composed of a hard core of diameter a which has 
a repulsive shoulder of width b — a at depth Ae, and an attractive well of width c — b and 
depth e. The form of the function is thus 

oo < r < a 
— Ae a < r < b 

(1) 

— e b < r < c 
r > c 



u(r) 



All of the results reported here for the discrete potential are for a = 1, b = a/2, c = \/3, 
e = 2, and A = 0.5. 

In the case of water, one can attribute the larger distance r = b to hydrogen bonding, 
for which the system acquires a low energy and expands at the same time. The inner 
distance r = a on the other hand corresponds to a non-hydrogen bonded energy state. 
Recent studies have proposed this form of potential as the interaction between clusters of 



strongly bonded pentamers of water This type of interaction is expected to reproduce 
the density anomaly. The reason is that at low pressures and at low temperatures, nearest 
neighbor pairs sit in the outer well which has a lower energy. By increasing T, in order to 
gain more entropy, the system explores a larger portion of the configurational space, which 



is not probed at lower temperatures. This includes penetrating the "soft" core (a < r < b), 
which on average can lead to anomalous contraction upon heating. 



B. Smooth Potential 



While the discrete potential u(r) is appropriate for deriving the closed form of equation 
of state in Id |j, |1^, [HJ , for simulations it is not necessarily the most appropriate. As we 
will show in the next section, the smooth version of the potential u'(r) requires a different 
method of simulation from that of u(r). The potential u'(r) we use is obtained by adding a 
Gaussian well to the Lennard- Jones potential and has the form 



u'(r) = 4e' 



X'e' exp 



a a 



(2) 



for r < r c and vanishes for r > r c . We use e' = 1.0, A' = 1.7, w = 5.0, tq = 1.5, a — 1, n — 2 
in order to mimic the shape of the discrete potential, as shown in Fig. |l|. 



III. MOLECULAR DYNAMICS SIMULATION METHODS 

The method of simulation in both the discrete and the smooth cases is the molecular 
dynamics (MD) method. Our simulations are performed in 2c? with periodic boundary 
conditions. The overall qualitative results of the simulations for the discrete and the smooth 
potential are similar, while the quantitative results differ. In what follows, we explain in 
more detail the MD method used in each case. 



A. Constant- Volume Simulation of the Discrete Potential 



For the discrete potential [Eq. ([[])], we use the collision table technique jnj for N = 896 
disks. To each disk we assign a radius a/2. We define the density p to be the ratio of 
the total area of all the disks to the area of the box. Energies are measured in units of e, 
temperature is calculated in units of energy divided by the Boltzmann constant, e/ks, and 
the mass of the particle is set at m — 1. 



The average pressure is calculated using the virial equation for step potentials fll4 



V 



(3) 



where K = J2i m ^i/^ is the total kinetic energy, N is the number of particles, J2'ij is the 
sum over the particle pairs undergoing a collision in the time interval St, Vi and v[ are 
the velocities of the particle % before and after a collision, r*j and fj are the positions of the 
particles i and j undergoing a collision at the start of St. 

We simulate state points along constant-volume paths. For thermalization we use the 
Berendsen method of rescaling the kinetic energy [14]. We thermalize the system for 10 5 
time units, which corresponds to ~ 10 6 collisions per particle, and then acquire data for 10 6 
time units corresponding to ~ 10 7 collisions per particle. 

B. Constant- Volume Simulation of the Smooth Potential 



For the smooth version of the potential [Eq. (|2])], we use the velocity Verlet integrator 



method |TJ] for a system of iV = 2500 discs. We record the results in reduced units in which 
a, e', m, and ks are all unity. We choose r c = 2.5, and the length of each MD time step 
St = 0.01. We assign to each particle a radius 2 1//6 , which corresponds to the minimum of 
the Lennard- Jones potential a, and define the density p to be the ratio of the total area of 
all the disks to the area of the box. 

In order to achieve a preset temperature we use the Berendsen method of rescaling the 



velocities [I4"|, resulting in the time dependence 



T{t) = T oc + (T(0) - T 00 )e- 



-t/r 



(4) 



where r is a preset time constant JT^[. Typical values of r are around 10 4 5t. 

We first thermalize the system for a time At w lOr, and we ensure that equilibrium 
is attained by monitoring the time dependence of observables like T, P and the potential 
energy U. Then we acquire data, running the system for an additional period of time at 
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constant NVE conditions (micro-canonical ensemble). We calculate P and T and we consider 
the system to be in equilibrium only when the fluctuations of these quantities are less than 
1% of their average values. The acquisition time is chosen to be more than the time it takes 
for the system to equilibrate and is typically (5 x 10 4 )<5t to (2 x 10 5 )<5t. 

The average simulation speed on Boston University's SGI Origin 2000 supercomputer 
was approximately 10/zs per particle update. Each of the state points we study requires 
between 8 and 16 hours on one processor, and thus over 1000 hours total computational 
time was utilized. 

C. Constant-Pressure Simulation of the Smooth Potential 

In order to check that our results are not biased by problems like phase coexistence that 
can affect constant-volume simulations, we also perform constant-pressure simulations in 
the case of the smooth potential. Constant-pressure simulations allow us to determine more 
accurately the locations of the freezing line and the density maximum. We use the feedback 
method proposed by Broughton |15|, where the dimensions of the box and the positions of 
the particles are rescaled to obtain the desired pressure P. The amount of rescaling depends 
on the difference between the present pressure P(t) and the desired pressure P. 

Using the Broughton method and the Berendsen method, we gradually drive the system 
to the desired P and T, while simulating under readjusting V and E conditions. We choose 
pressure and temperature tolerances 5P and 5T less than 1% of the desired P and T. Once 
P(t) and T{t) reach values within the range P ± 5P and T ± 5T, we stop thermalization 
and pressurization. If the system stays within these limits for an interval of time of the 
same order of time needed to reach the desired P and T, we conclude that the system 
has equilibrated, turn off the thermalization and pressurization, and start collecting data. 
During this collection period, we monitor temperature and pressure to check that their 
average values coincide with the desired ones within an accuracy of 1%. For our results, 
the time At needed to reach equilibrium is usually of the order of 500,000 steps 5t, so 
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At ~ 5000 in Lennard- Jones units. Data are acquired over a period of lOAt. We test our 
code by simulating a Lennard- Jones system of 2304 disks and comparing the results with 
the previously-known phase diagram of a 2d Lennard- Jones system. 



IV. DENSITY ANOMALY 

The temperature of density maximum (T^/d) line is the border of the region in the P — T 
plane where the liquid expands upon cooling. Fig. shows a set of different isobars for the 
smooth potential. The Tud line corresponds to the set of maxima along those isobars. 

In the case of constant- volume simulations, the Tmcl line corresponds to the minima along 
constant-volume paths (P vs. T graphs of Fig. 3) since for any thermodynamic quantity X 
constant-volume paths 

(dX\ (dX\ {dX/dP) T {dV/dT) p 

\dT) v \dT) p (dV/dP) T ■ { > 



By substituting X = P we find 



'dP\ 

dT) y K T 



) (6) 



where 



dT J 



is the thermal expansion coefficient and 



is the isothermal compressibility. Taking a derivative of Eq. (||), using ap = at the T md 
line, we find 



' d 2 p \ v- 1 (d 2 v/dr 2 ) p 

dT 2 )„~ K T 



(9) 

i \ - 1 ■ 

Equations (|5]) and (Q) show that since is always positive and finite for systems in 
equilibrium not at a critical point, a minimum along the isochore is equivalent to a minimum 



along an isobar, which is the density maximum point Tud- Figure |3] shows the isochores for 
the smooth and the discrete potentials. 

To confirm that we are investigating the liquid state part of the phase diagram, we 
introduce a criterion to distinguish the liquid state from a frozen state. We determine the 
freezing line as the location of points where isochores overlap. In this way we establish an 
approximate location for the freezing line. Crossing this line from the liquid side, we find 
a sharp decrease in D coinciding with the appearance of slowly-decaying peaks in g(r) as 
a function of r, which signals the build up of long-range correlations (Fig. f|), which is a 
characteristic of 2D solids. 

We confirm the above criterion adopted to locate the freezing lines by using isobaric 
simulations for the smooth potential. Indeed, they show a sharp change of density, in 
correspondence with the estimated freezing line at high pressures (Fig. [|, lower panel). The 
presence of a hysteresis loop (Fig. ||) suggests that the liquid-solid transition is first-order; 
however, by lowering the pressure the loops become less and less pronounced and eventually 
disappear, leaving the possibility of an hexatic second-order transition |TB| . 

For a few state points near the freezing line, we have checked our results by simulating 
N = 2500 particles in rectangular boxes with aspect ratio Vo/2 (L x = \^3/2L y , with 
L x x L y = Lq) which accommodate triangular lattices perfectly. This eliminates any possible 
artificial hindrance in crystallization that may arise from the asymmetry imposed by the 
shape of the square box. 

For water, the locus of the Tm& line i n the P-T phase diagram is of special interest to 
distinguish between different scenarios proposed to explain its anomalies [pTTl — jTSfl - In Fig. [7| 
we see that the Tud line changes from negative slope at high pressures to positive slope 
at low pressures. This change in slope is similar to what is observed in simulating model 



potentials of water like SPC/E or ST2 |T7|,pO[. These results suggest that the change in 
slope can be a general phenomenon stemming from the general core-softened form of the 
interaction in the simulation. 
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V. STRUCTURES IN THE LIQUID AND SOLID PHASES 



Figure |5] (upper panel) shows the different phases of the system. In our simulations 
we see that the Tud line is located in the region of pressures where the freezing line is 
negatively sloped, as in water. A density anomaly and a negatively-sloped melting line are 
often associated HHU- This has proven to be the case for substances like water (Fig. H) 



and tellurium M and for computer models 12,22]. This association is plausible since the 



isobaric thermal expansion coefficient otp is related to the cross fluctuations in volume and 
entropy as 

a P = (3(5V5S). (10) 

Approaching a freezing line, we expect local density fluctuations to have structures similar 
to the neighboring solid as they are going to trigger the liquid-solid transition. On the other 
hand, the Clausius-Clapeyron relation for the slope of the freezing line 

^ = ^ (11) 

implies that, if the freezing line is negatively sloped, the solid, which has a lower entropy 
than the liquid, will have a higher specific volume. Therefore, if the fluctuations in the liquid 
are "solid-like," ap [Eq. (|I"0"D1 will turn out to be negative. 

To distinguish different local structures in the liquid, we plot the radial distribution 
function g(r) for different pressures and temperatures (Fig. ^). As expected, at low pressures 
cooling expels particles from the core, while increasing pressure at fixed temperature has 
the opposite effect. 

Since our system is two-dimensional, we can use visual inspection to develop an intuitive 
picture of the possible local structures (Fig. |, upper panel). If the fluctuations in the liquid 
are "solid-like," near the freezing line we expect to see local structures that resemble the 
structure of the nearby solid. 

We find that at low P and T, the system is frozen with a hexagonal structure (Fig. 
lower left panel). A "snapshot" of the system along the same isobaric line (Fig. || lower right 
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panel) shows clearly that local patches with hexagonal order are present in the liquid phase 
near the freezing line. We will refer to this structure as the "open structure." Similarly, at 
high pressures the local patches in the liquid phase near the freezing line (Fig. [^, upper right 
panel) resemble the structure of the system when it is frozen at low T and high P (Fig. |9|, 
upper left panel). We will refer to this as the "dense structure." For the open structures, 
each particle has six neighbors sitting in the deepest well, and the softened-core behaves as 
the effective core for the particles. The dense structure is the next energetically favorable 
local arrangement, with four neighbors in the external well and four in the softened core, 
for which the effective core is the hard core. 



VI. DIFFUSION ANOMALY 

We next study the diffusion anomaly, which is another surprising feature of water. While 
for most materials diffusivity decreases with pressure, liquid water has an opposite behavior 
in a large region of the phase diagram || (Fig. ^[). The pressures where the system has a 
maximum diffusivity along isotherms define the line of the pressure of maximum diffusivity, 
Pmd- 

We observe that our core-softened potential reproduces this anomaly. We first mea- 
sure the mean square displacement (Ar 2 (t)) = ([r(t + t ) — r(t )} 2 ) and then the diffusion 
coefficient using the relation 

D = 1 lim <^!M. (12) 

We measure (Ar 2 (£)) by averaging over the starting time to in Eq. ([12]). We find that there 
is a region of the phase diagram in which D increases upon increasing P (Fig |10D . 

In order to understand the diffusion anomaly we first note that, for normal liquids, D 
decreases with P because upon increasing P the density increases and molecules are more 
constrained. As a result the viscosity increases and hence D decreases. In the case of 
water the anomaly can be related to the fact that increasing pressure (and hence density) 
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breaks hydrogen bonds, which in turn increases the mobility of the molecules. We present a 
more general explanation which can equally apply to our radially symmetric core-softened 
interaction which does not possess any directional bonds similar to hydrogen bonds. The 
low energy inter-particle state at r w b plays the role of non-directional bond. Note that 
D is proportional to the mean free path of particles, which increases with the free volume 
per particle Vf ree = v — v ex , where v ex is the excluded volume per particle resulting from the 
effective hard core. At low temperatures, v ex for the dense structure is proportional to the 
area a 2 of the hard core, while for the open structure it is proportional to the area b 2 of 
the soft core. Increasing P decreases v , which is the main effect in normal liquids. For the 
core-softened liquid, on the other hand, increasing P can also decrease v cx by transforming 
some of local open structures to dense structures. Since both Av and Av ex decrease with 
P and since Avf ree = Av — Av ex , the effect of P on D depends on whether Av or Aw ex 
dominates. The anomalous increase in D along the isotherms near the freezing line is a 
sign of the dominance of the Ai> ex term. Thus the anomaly in D must disappear near a 
certain pressure above which the average distance between particles corresponds to the dense 
structure, and as a result the contribution of the open structure to t> ex is negligible. 

We verified this in our simulation by observing a correspondence between the disappear- 
ance of the diffusion anomaly and the disappearance of the peak in g(r), corresponding to 
the open structure that is observed in real water. 



VII. ISOTHERMAL COMPRESSIBILITY 

In order to investigate the anomaly in isothermal compressibility Kt, we calculate Kj- 
at each state point using the data in Fig. 3. In the smooth potential case we verify these 
results using 

K T = \im^^- (13) 



as an alternative method [53], where n is the density of the system and S(q) is the structure 



function and is related to the pair correlation function via 
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S(q) = l + nj e iC l- x ^(x)rfx. 



(14) 



We first calculate g(r) for each state point averaging over all thermalized configurations. 
We then perform numerical integration using Eq. ( |i~4"D to find S(q), and finally we extrapolate 
S(q) to q = and substitute in Eq. (|T3| ) to find Kt- We show an example of the graphs for 
g(r) and the resulting S(q) (normalized by the extra factors in Eq. [IS]) in Fig. O. We fit 
the low q tail of the curve to a line to find the limiting value in Eq. (|i"3[). 



We graph Kt along isobars, as shown in Fig. [12]. For large T, the ify decreases upon 
increasing P. For small T the behavior is the opposite and the compressibility anomaly 
occurs. As seen for all isobars shown in Fig. [12] (except a low pressure one), Kt increases 
by lowering T. 

VIII. PHASE DIAGRAM 

In water the T^d line is negatively sloped for positive pressures. For several models 
that mimic water behavior, it is found that the Tud line has a reentrant shape, changing 
slope at low or negative pressures |3J . In our simulations, we find such a reentrant Tud line; 
the change of slope of the Tud happens at positive pressures in the smooth version and at 
negative pressures in the discrete case (Fig. |7|). 

Moreover, we have graphed the location of the minimum Kt point along each isobar; 
the locus of these points is called the minimum compressibility line [Kt ■ )• Sastry and 



coworkers | 18[] , from basic thermodynamic arguments, show that: (i) the Kt . line intersects 
the TMd line at its infinite slope point, and (ii) the compressibility must increase upon cooling 
in the region to the left of a negatively sloped Tj^d line. Our results are in agreement with 
both of these statements (Fig. 0). 

Theories relating D to the entropy pi] predict that the anomalous behavior (dD / 8P)t > 
is related to an anomaly in the entropy (dS/dP)T > 0. Due to the Maxwell relation 
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(dS/dP)T = —(dV/dT)p, whenever there is a density anomaly, an entropy anomaly occurs, 
and the value of entropy along isotherms reaches a maximum at the line. 

In Fig. |7| we also show the Pmd line where D reaches its maximum with pressure. Notice 
that, for the continuous potential, the maximum in D shifts to higher P with increasing T. 
This trend is also observed in the SPC/E model of water [35] but is in contrast with the 
behavior of real water (Fig. ||). 



IX. SUMMARY 

We find that core-softened potentials reproduce the qualitative behavior of water in 
many respects; in particular, the liquid phase of core-softened potentials can show both 
thermodynamic anomalies and dynamic anomalies. Moreover, as in real water, the freezing 
line changes slope from a positive value at high pressures to a negative value at low pressures 
in the P-T phase diagram and more than one solid phase is present. The polymorphism 
of the solid phase and the anomalies in the liquid phase can be related to the possibility 
of different local structures due to the shape of the potential. The phase diagrams of the 
discrete and the smooth versions of the core softened potential are similar to that of real 
water, but the T Md line is shifted into liquid phase and the K Tmin line has a positive slope. 
However, only for the discrete potential we find a Pmd line with a negative slope. 
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FIGURES 
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FIG. 1. Discrete (broken line) and smooth (solid line) forms of the core-softened potential 
studied here. Length parameters a, b, c and energy parameters e, A are shown. 
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FIG. 2. Isobars for the smooth potential 
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FIG. 3. Isochores for the smooth potential (left panel) and the discrete potential (right panel). 




FIG. 4. Radial distribution function at high and low temperatures, along the P=0.48 isobar 
for the smooth potential. Notice how, by lowering T, long range correlations develop (g(r) = 1 if 
particles at distance r are uncorrelated) and more particles are expelled from inside the soft core 
r ~ 1.1 into the attractive well r ~ 1.5 (inset). As the average interparticle distance is growing 
upon cooling, the system is expanding and there is hence a density anomaly. 
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FIG. 5. Phases for the core-softened model (smooth potential). The upper panel shows a 
snapshot of the liquid phase and snapshots of different types of crystals for the solid phase at low 
pressures where the freezing line is negatively sloped, and at high pressures where the freezing line 
is positively sloped. Lower panel shows the density jumps along isobars; note that the low pressure 
isobar shows a density anomaly before jumping to a low density solid. 
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FIG. 6. Hysteresis loop near the freezing line for a high pressure isobar (smooth potential). 
The continuous line is obtained upon heating, the dashed line upon cooling. 
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FIG. 8. Sketch of the phase diagram of water. The portion of the T^d line that is to the left 
of the melting line corresponds to experiments in the supercooled region of water. Notice that the 
presence of a density anomaly in the region of the negatively sloped melting line can occur in the 
metastable phase of the liquid. Data are obtained from Refs. 0. 
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FIG. 9. Snapshots of the system (discrete potential) in the solid phase at high pressure (upper 
left) and low pressure (lower left), and in the liquid phase at high pressure (upper right) and low 
pressure (lower right). Moving along an isobar, patches of local order similar to the low temperature 
solid develop. At high pressure, the average distance between particles (which is the radius of the 
disk) is of the order of the hard core, while in the low pressure solid, the distance between particles 
is larger, of the order of the softened core. 
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FIG. 10. Diffusion coefficient D in the liquid phase for the continuous potential (left panel), 
the discrete potential (right panel) along various isotherms. Lines are intended as a guide for the 
eye. Notice the anomalous sections of the graph, where (dD j '8P)t > 0. 
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FIG. 11. Averaged pair distribution function for the smooth potential at P = 1.0, T = 0.7 
(upper panel). The structure function (lower panel), multiplied by the factor 1/nksT, derived by 
integrating g(r), where ks is the Boltzmann constant. The q — > limit of this function gives Kt, 
which is around 0.1 in this case. 
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